An Example of How To Design, Size, and Test a Vehicle Thermal Management System

Keywords: thermal system sizing, thermal system design, thermal system digital twin, BEV or HEV system-level design, energy management studies
Target audience: Automotive engineers that are involved in system-level design and are not familiar with refrigeration or temperature-controlled system design (for instance, engineers with a mechanical or electrical background)
Table of Contents

System Overview

In this example, we expand the previously designed simple refrigeration loop system (for details, refer to designHVACsystem.mlx Live Script) with additional components, for instance the vehicle cabin and a chiller that interfaces with a water-based cooling circuit, thus creating a more realistic model of a vehicle thermal management system. Here too, we use energy flow based system-level heat exchangers, and not geometry-based heat exchangers. Usually, the OEM system-level engineer has little knowledge regarding detailed heat exchanger design and the OEMs will source them from their suppliers based on power requirements.
The thermal management system uses a standard vapor compression cycle (the refrigeration loop) that is used for direct cabin cooling and an indirect battery cooling via a secondary water-glycol loop. Superheated R410a refrigerant is compressed across a parameterized compressor model and then condensed via heat exchange with ambient air. The refrigerant flow splits and can expand across two expansion valves. The first valve (TXV1) leads to the cabin heat exchanger while the second valve (TXV2) leads to another heat exchanger which acts as a chiller and cools a secondary water-glycol flow loop that conditions the electric components (like battery) .
The contribution's target audience is an automotive engineer that is involved in system-level design and is not familiar with refrigeration or temperature-controlled system design.
Hopefully this example will help them understand how to properly size and design a system-level thermal system model for a BEV, HEV, or a conventional vehicle.
After going through the example, the engineer will be able to properly size, design, and test the thermal system of a vehicle using a virtual testbench (aka digital twin).
The goal of the model is to be used as a virtual testbench that can be used to determine the size of the system components (and then, for instance, using supplier parts catalogues to select appropriate components).
This model can be used to properly size the thermal management system various components. The main interest is in plant modeling, understanding the heat energy rates in the various heat exchangers and the power consumed by compressor, fans, and fluid pumps.
The focus is not on developing control algorithms, and for this reason the control algorithm is minimal, the only feedback control used is a bang-bang control for cabin temperature.

System Components

Refrigeration System

NOTE: To understand how the refrigeration loop parameters used below were defined, refer to designHVACsystem.mlx Live Script for a detailed explanation.
% Set nominal cooling temperature of the refrigeration system
T_set = 22; % [degC] Thermostat set point
 
% Cooling capacity of the refrigeration system (this parameter defines the
% size of the system)
CoolingPower_kW = 3; % [kW]
 
% Thermodynamic design parameters
evap_inlet_enthalpy = 267; % [kJ/kg] (from p-H diagram)
evap_outlet_enthalpy = 430; % [kJ/kg] (from p-H diagram)
 
evap_inlet_p = 0.934; % [MPa]
evap_superheat_temp = 5; % [deltaK]
 
% Set the evaporating temperature, or the saturation temperature in the
% evaporator, to be lower than the desired cabin temperature to enable
% heat transfer from the indoor air to the refrigerant.
nominal_evap_temp = T_set -15; % [degC]
 
% Calculate the refrigerant mass flow rate using enthalpy difference
refrigerant_massflowrate = CoolingPower_kW/(evap_outlet_enthalpy-evap_inlet_enthalpy); % [kg/s]
 
% To calculate air massflow rate, divide the cooling capacity by the heat
% capacity of air at constant pressure, and divide that result by the
% desired temperature drop across the evaporator
Cp_air = 1.005; % [kg/kJ/K] Cp heat for air
evap_temperature_drop = 10; % [degC], temperature drop across the evaporator
evap_air_massflowrate = CoolingPower_kW/Cp_air/evap_temperature_drop; % [kg/s]
 
% Define refrigerant pipe and air duct diameters
tube_D = 0.01; % [m] Refrigerant tube diameter
duct_W = 0.2; % [m] Air duct width
 
% Set the condensing temperature, or the saturation temperature in the
% condenser, to be higher than the outdoor temperature to enable heat
% transfer from the refrigerant to the outdoor environment.
nominal_condens_temp = 45; % [degC]
% Set a subcooling temperature of 5 deg
condenser_subcooling_temp = 5; % [deltaK]
condenser_inlet_temp = 65; % [degC]
condenser_inlet_enthalpy = 457; % [kJ/kg (from p-H diagram)
condenser_inlet_p = 2.734; % [MPa] (from p-H diagram)
 
% Mass flow rate for the condenser fan (design parameter)
condenser_air_massflowrate = 0.5; % [kg/s]
 
% Compressor nominal speed
RPM_Compressor_nominal = 1000; % [rpm]

Vehicle Cabin (HVAC volume)

The cabin is modelled as a masked subsystem, as below.
The subsystem contains a thermal network of the cabin and a moist volume.
The parameters are defines as below:
cabin_L = 2; % [m] cabin length
cabin_W = 2; % [m] cabin width
cabin_H = 1.5; % [m] cabin height
door_or_window_area = 0.5; % [m^2] Door or window area
T_cabin_init = 30; % [degC] Initial cabin temperature
RH_cabin_init = 0.5; % [-] Initial cabin relative humidity
 
The cabin is heated by the solar irradiation and exchanges heat with the environment and with the cooling circuit.
Heat convection with the environment depends on the vehicle speed.
Cabin temperature is controlled using the HVAC evaporator and the blower fan.
The following figure shows the vehicle cabin model. The subsystem contains a moist air volume, thermal network of the cabin, and an occupant heat model. The moist air volume models the cabin air.
The cabin's thermal network is a detailed model with the parameters defined in the block mask. Cabin glass, doors, roof, and the occupants thermal model (human model), are modelled using a Simscape thermal network and a moist air source.
Cabin's thermal network model is defined as below:

Condenser, Evaporator, and Coolant Radiator fans

Condenser and evaporator fans are modelled with a variant subsystem, having two levels of fidelity:
  1. Simple fan model, using a volumetric mass flow rate block
  2. Detailed fan model, using the Fan (MA) block available from Simscape Fluids product.
Fan parameters are defines as below:
air_density = 1.225; % [kg/m^3] at standard conditions
VolFlowrate_CondenserFan = condenser_air_massflowrate/air_density;% [m^3/s]
RPM_CondFan = 3000; % [rpm]
VolFlowRate_EvapFan = evap_air_massflowrate/air_density;% [m^3/s]
RPM_EvapFan = 3000; % [rpm]
VolFlowRate_CoolantRadiatorFan = 0.5;% [m^3/s]
As for the refrigeration compressor, since it is an important element that can maintain a stable operation of the refrigeration cycle, in this case we will use the detailed model variant. One thing to note here, during operation the refrigerant flow rate is lower bounded to a value slightly bigger than zero. That is why we limit the compressor's minimum speed at 10 rpm, refer to below figure. Maintaining a very small non-zero flow rate will enhance simulation robustness. Additionally, since the compressor control is for the moment just a bang-bang control, we use a low-pass filter with a time constant of 2 s to smooth out the compressor on/off command value. This will help the numerical solver take larger steps during simulation. Such input smoothing is used for other system actuators as well, enabling faster model initialization and simulation.

Chiller for the Cooolant Network

The chiller is a heat-exchanger that functions as an evaporator, and which couples the refrigerant loop with the water-based coolant circuit.
With the exception of fluid flow mass flow rate, the parameters are the same as the ones used for the HVAC evaporator, if needed they can be adjusted to resize the system.
The thermal liquid nominal mass flow rate is calculated based on the coolant pump displacement, see further down for calculation.
Just like the evaporator, the chiller refrigerant flow is also controlled with a Thermostatic Expansion Valve (TXV) whose parameters are defined as below:

Coolant Circuit Components

This subsystem models a water-based cooling circuit, that can be used to control the temperature of various electrical devices, such as batteries, DC/DC converters, inverters, motors, etc.
For the moment, the purpose is system-level estimation of the cooling power needed, so the coolant circuit is a simple model. To increase the level of fidelity, the coolant circuit can be interfaced with a more detailed system that includes a BEV or HEV electric components or an ICE (internal combustion engine). This part is out-of-scope for the present example, and the reader can consult MathWorks® documentation, for example Electric Vehicle Thermal Management.
To size the thermal system, for the moment we simply assume a heat flow as input. Typically, this is a known quantity that can be easily calculated summing up the various components heat losses.
The heat flow is exchanged with the coolant circuit using a pipe element with parameters defined below:
coolant_pipe_D = 0.02; % [m] coolant pipe diameter
coolant_circuit_thermal_mass = 5; % [kg]
coolant_circuit_length = 5; % [m]
The coolant circuit contains a fan radiator that is used to complement the chiller.
Radiator parameters are defined as below:
Radiator fan is modelled as a simple volumetric flow rate source, and radiator_air_mass_flow_rate is calculated below.
radiator_cooling_power = 3; % [kW] radiator cooling power
air_density = 1.225; % [kg/m^3]
radiator_coolant_mass_flow_rate = 1; % [kg/s]
radiator_air_mass_flow_rate = air_density*VolFlowRate_CoolantRadiatorFan; % [kg/s]
The coolant pump and tank parameters are defined as below:
pump_displacement = 0.02; % [l/rev] coolant pump displacement
pump_speed_max = 3000; % [rpm] coolant pump rpm
coolant_density = 1092; % [kg/m^3] Density of Ethylene Glycol in Solution for Mass Fraction of 0.5
pump_nominal_massflowrate = pump_displacement*1e-3*(pump_speed_max/60)*coolant_density; % [kg/s] used to parametrize the chiller's thermal fluid nominal mass flow rate
The Fixed-Displacement Pump is parametrized based on its displacement.
Coolant tank is modelled using Tank (G-TL) block, which is a pressurized tank with variable gas and thermal liquid volumes.
For simplicity, coolant tank is assumed not to exchange heat with the surroundings (if needed, one can add a simple thermal network to enable heat exchange). Heat exchange with the surroundings can play an important role, since it can be rather large due to high heat convection occurring at higher vehicle speeds.
coolant_tank_volume = 5; % [l] coolant tank capacity
coolant_tank_area = 0.11^2; % [m^2] coolant tank area

Initial Parameters

% Environment
T_env = 30; % [degC] External environment temperature
RH_env = 0.5; % External environment relative humidity
 
moist_air_p_init = 0.101325; % [MPa]
moist_air_T_init = T_env; % [degC]
moist_air_RH_init = RH_env;
 
coolant_p_init = 0.101325; % [MPa] initial ccolant pressure
coolant_T_init = 40; % [degC] initial coolant temperature
 
% Initial refrigerant condition is set up automatically in the block mask
% using nominal values, if not uncomment and use below values
% refrigerant_T_init = 20; % [degC] initial
% refrigerant temperature
% refrigerant_p_init = 1.0; % [MPa] initial refrigerant pressure
% refrigerant_alpha_init = 0.76;% initial refrigerant vapor quality
We are now ready to run the model.

Simulation

We assume a simple test scenario described below:
modelName = "VehicleThermalSystemSizing";
open_system(modelName);
First, let's investigate the performance of the model using a variable-step solver.

Simulation with a variable-step solver

We will first use a variable-step solver for the simulation: daessc.
The daessc variable-step Simulink® solver provides algorithms specifically designed to simulate differential algebraic equations (DAEs) arising from modeling physical systems.
set_param(modelName,SolverName="daessc",MaxStep="1");
% set heat flow rate to the coolant to be 0.5 kW
set_param(modelName+"/Heat flow rate to coolant [kW]", Value = "0.5");
To check the Simscape model initialization (sometimes there are issues with initialization), first we will enable Simscape data logging (in case we need to debug the Simscape model), and collect simulation data for the first 20 seconds. We also comment on the p-H diagram, to increase the simulation speed.
set_param(modelName,...
SimscapeUseOperatingPoints = "off",...
SimscapeLogType = "all", ...
SimscapeLogName ="simlog_refrigeration",...
SimscapeLogOpenViewer = "off");
 
set_param(modelName+"/P-H Diagram", Commented = "on" );
 
% Simulate the model for 20 seconds.
out = sim(modelName, StopTime = "20");
The initialization time is not long and results looks good but the first few seconds of the simulation are slow, indicating that the variable-step solver needs to take very small steps.
This may cause issues when trying to solve the model with a fixed-step solver, we will discuss this issue later.
Next, simulate the model for 300 seconds.
out = sim(modelName, StopTime = "300");
Once again, the results look good. The thermostat is keeping cabin temperature around set value of 22 °C and the coolant temperature goes down, which show the chiller and radiator are working. The energy flow rates look realistic. The simulation is stable even if we switch on and off the compressor.
When the compressor is on, condenser heat flow rate is around 3 kW, and the cabin temperature reaches 22 ℃ in approximately 1 min, as specified in the design requirements of the refrigeration loop.
The cabin temperature becoming lower than the outside temperature, the heat flow rate exchanged with the exterior reaches almost -500 W, minus sign indicating that the cabin is heated by the hotter environment.
Coolant temperature at steady-state is around 31 ℃.
Let's see how fast the model is running, by using out.SimulationMetadata.TimingInfo.
Notice that the simulation speed is relatively high, a few times higher than wall clock time.
Next, let's investigate the performance of the model when using a fixed-step solver.

Simulation with a fixed-step solver

Such solvers must be used when deploying the model to a real-time target, for instance a Hardware-In-the-Loop (HIL) platform.
We use the ode1be solver with a step time set to 0.2 seconds. Step time is rather big, but we need this for a higher simulation speed.
Hopefully the simulation does not become unstable with such a large step.
set_param(modelName,SolverName="ode1be", FixedStep = "0.2");
To check the Simscape model initialization, first we will enable Simscape data logging (for debugging), and collect simulation data for the first 20 seconds.
set_param(modelName,...
SimscapeUseOperatingPoints = "off",...
SimscapeLogType = "all", ...
SimscapeLogName = "simlog_refrigeration",...
SimscapeLogOpenViewer = "off");
 
% Simulate the model for 20 seconds.
out = sim(modelName, StopTime = "20");
This time the simulation speed for the first few time steps or so is very low! It takes almost a minute to execute the first steps.
This was indicated as a potential issue when we have run with the variable step solver.
Notice that there are some initial transients in the first few seconds when using the fixed-step soolver, indicating possible instability or numerical issues.
We may try to look at the initial parameters and see if there is any problem, or try to lower the fixed-step time - but this will slow down the simulation speed.
We will pursue a different approach.
Another solution is to use a Simulink.op.ModelOperatingPoint object and use it to initialize the model before another simulation run. Below is an example how to do this.

Get and Set initial state of the model

% Create an operating point from logged simulation data at 1 seconds after the start of simulation:
in1 = Simulink.SimulationInput(modelName);
in1 = setModelParameter(in1,StartTime="0",StopTime="1",SaveFinalState="on",SaveOperatingPoint="on");
 
out1 = sim(in1);
 
% Enable model initialization using final state out1.xFinal, and run the
% simulation from 1 s to 300 s
in2 = Simulink.SimulationInput(modelName);
in2 = setModelParameter(in2,StopTime="300");
in2 = setInitialState(in2,out1.xFinal);
 
% Run the simulation from 1 s to 300 s
% Note that the initial condition is inherited from the operating point
% used (for instance the initial coolant temperature is not 40 degrees but
% slightly lower).
out2 = sim(in2);
The simulation starts immediately and the model runs pretty fast!
Notice the simulation starts at 1 s and not 0 s.
Let's check the execution time.
Results indicate that the model runs at almost 3x faster than clock time, so with some care it may run on real-time on a HIL system, for instance.
Next, let's try to see what happens if the heat flow rate from the electrical components (battery, motors, etc.) is changed from 0.5 kW to 2 kW.
in2 = setBlockParameter(in2,modelName+"/Heat flow rate to coolant [kW]",Value="2");
 
% Run the simulation from 1 s to 300 s
% Note that the initial condition is inherited from the operating point
% used (for instance the initial coolant temperature is not 40 degrees but
% slightly lower).
out2 = sim(in2);
 
This time the coolant temperature is higher than before, as expected due to the higher rate of input heat flow. However, the chiller and radiator fan maintain a coolant temperature of about 34 ℃, and the refrigeration system can still keep the cabin cool.

Conclusions

What have we shown in this demo? We have shown how to build a system-level digital twin for a vehicle thermal management system.
  1. We have started by designing a refrigeration loop using thermodynamic principles
  2. We expanded the model to handle an additional water-based cooling network
  3. We investigated simulation speed and how to solve model initialization issues
Based on the energy flow considerations, you may need to adjust the cooling power of the thermal system to match your downstream coolant-based circuit, to increase the cooling rate for faster cabin cooling, or to satisfy cost vs. performance design constraints. This is where such system-level models come in handy, enabling detailed yet fast what-if scenarios!

Dependencies

MathWorks products needed to run the simulation are listed below:
  1. MATLAB® version R2023b
  2. Simulink®
  3. Simscape™
  4. Simscape™ Fluids™

License

The license used in this contribution is the XSLA license. Refer to LICENSE for more information.
Copyright 2024 The MathWorks, Inc.